Podsumowanie analizy

Przedmiotem analizy był zbiór zawierający pomiary trzyletnich śledzi oraz stanu środowiska, w którym żyją. Dane były zbierane na przestrzeni ostatnich 60 lat. Celem analizy było odkrycie przyczyny spadku długości śledzi na przestrzeni lat.

W ramach analizy dokonano podsumowania rozkładów zebranych atrybutów. Policzono także korelacje między zmiennymi. Okazało się, że najbardziej skorelowane są atrybuty reprezentujące liczbę różnych gatunków planktonu. Najmocniej skorelowana z długością śledzi jest natomiast temperatura przy powierzchni wody. W połączeniu z analizą nauczonego regresora prawdopodną wydaje się hipoteza, że atrybut ten w największym stopniu przyczynił się do obserwowanego spadku liczby śledzi.

Wykorzystane biblioteki

library(corrplot)
## corrplot 0.84 loaded
library(mice)
## Loading required package: lattice
## 
## Attaching package: 'mice'
## The following objects are masked from 'package:base':
## 
##     cbind, rbind
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(caret)
library(reshape2)
library(knitr)
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

Powtarzalność wyników

W celu otrzymywania powtarzalnych wyników ustalono stałą wartość ziarna, która zostaje wykorzystywana przez algorytmy korzystające z pseudolosowości:

seed <- 420
set.seed(seed)

Wczytywanie danych z pliku

Wczytywanie danych odbywa się z użyciem funkcji read.csv. Przekazanie znaku zapytania jako argumentu na.strings powoduje, że wszystkie brakujące wartości zostaną poprawnie oznaczone jako NA.

data <- read.csv("sledzie.csv", na.strings = "?")

Przetwarzanie brakujących danych

W celu uzupełnienia brakujących wartości wykorzystano bibliotekę mice. Skorzystano z metody cart, która stosuje drzewa regresyjne w celu predykcji brakujących wartości:

data <- complete(mice(data, meth='cart', seed = seed))
## 
##  iter imp variable
##   1   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   4  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   1   5  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   4  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   2   5  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   4  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   3   5  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   4  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   4   5  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   1  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   2  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   3  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   4  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst
##   5   5  cfin1  cfin2  chel1  chel2  lcop1  lcop2  sst

Podstawowe statystyki dotyczące zbioru

Rozmiar zbioru:

nrow(data)
## [1] 52582

Liczba atrybutów:

ncol(data)
## [1] 16

Przykładowe wartości (w tym przypadku pierwsze 6 pomiarów):

head(data)
##   X length   cfin1   cfin2   chel1    chel2   lcop1    lcop2  fbar   recr
## 1 0   23.0 0.02778 0.27785 2.46875 21.67333 2.54787 26.35881 0.356 482831
## 2 1   22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 3 2   25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 4 3   25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 5 4   24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
## 6 5   22.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831
##        cumf   totaln      sst      sal xmonth nao
## 1 0.3059879 267380.8 14.30693 35.51234      7 2.8
## 2 0.3059879 267380.8 14.30693 35.51234      7 2.8
## 3 0.3059879 267380.8 14.30693 35.51234      7 2.8
## 4 0.3059879 267380.8 14.30693 35.51234      7 2.8
## 5 0.3059879 267380.8 14.30693 35.51234      7 2.8
## 6 0.3059879 267380.8 14.30693 35.51234      7 2.8

Szczegółowa analiza wartości atrybutów

Z pomocą funkcji summary można podejrzeć minimum, maksimum, średnią oraz 1. i 3. kwartyl dla każdego atrybutu:

summary(data)
##        X             length         cfin1             cfin2        
##  Min.   :    0   Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000  
##  1st Qu.:13145   1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778  
##  Median :26291   Median :25.5   Median : 0.1111   Median : 0.7012  
##  Mean   :26291   Mean   :25.3   Mean   : 0.4460   Mean   : 2.0257  
##  3rd Qu.:39436   3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936  
##  Max.   :52581   Max.   :32.5   Max.   :37.6667   Max.   :19.3958  
##      chel1            chel2            lcop1              lcop2       
##  Min.   : 0.000   Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849  
##  1st Qu.: 2.469   1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808  
##  Median : 5.750   Median :21.673   Median :  7.0000   Median :24.859  
##  Mean   :10.002   Mean   :21.219   Mean   : 12.8084   Mean   :28.421  
##  3rd Qu.:11.500   3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232  
##  Max.   :75.000   Max.   :57.706   Max.   :115.5833   Max.   :68.736  
##       fbar             recr              cumf             totaln       
##  Min.   :0.0680   Min.   : 140515   Min.   :0.06833   Min.   : 144137  
##  1st Qu.:0.2270   1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068  
##  Median :0.3320   Median : 421391   Median :0.23191   Median : 539558  
##  Mean   :0.3304   Mean   : 520367   Mean   :0.22981   Mean   : 514973  
##  3rd Qu.:0.4560   3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351  
##  Max.   :0.8490   Max.   :1565890   Max.   :0.39801   Max.   :1015595  
##       sst             sal            xmonth            nao          
##  Min.   :12.77   Min.   :35.40   Min.   : 1.000   Min.   :-4.89000  
##  1st Qu.:13.60   1st Qu.:35.51   1st Qu.: 5.000   1st Qu.:-1.89000  
##  Median :13.86   Median :35.51   Median : 8.000   Median : 0.20000  
##  Mean   :13.87   Mean   :35.51   Mean   : 7.258   Mean   :-0.09236  
##  3rd Qu.:14.16   3rd Qu.:35.52   3rd Qu.: 9.000   3rd Qu.: 1.63000  
##  Max.   :14.73   Max.   :35.61   Max.   :12.000   Max.   : 5.08000

Wykresy pudełkowe dla każdego atrybutu (z pominięciem identyfikatora i miesiąca) prezentują się następująco:

data_m <- melt(select(data, -X, -xmonth))
## No id variables; using all as measure variables
ggplot(data = data_m, aes(x=variable, y=value)) + geom_boxplot() + facet_wrap( ~ variable, scales="free", ncol = 4)

Korelacja między zmiennymi

W celu odkrycia najbardziej skorelowanych zmiennych przetransformowano macierz korelacji w macierz trójkątną oraz wypisano wszystkie znaczące korelacje (o współczynniku korelacji większym niż 0.5):

cor <- cor(data, use = "complete.obs")
corUpper <- cor
corUpper[lower.tri(cor, diag = TRUE)] <- 0
corHighest <- subset(as.data.frame(as.table(corUpper)), abs(Freq) > 0.5)
corHighest <- corHighest[order(-abs(corHighest$Freq)),]
colnames(corHighest) <- c("Atrybut 1", "Atrybut 2", "Współczynnik korelacji")
corHighest
##     Atrybut 1 Atrybut 2 Wspólczynnik korelacji
## 101     chel1     lcop1              0.9562203
## 118     chel2     lcop2              0.8861733
## 169      fbar      cumf              0.8165530
## 187      cumf    totaln             -0.7077921
## 116     cfin2     lcop2              0.6532125
## 247     lcop1       nao             -0.5503011
## 253       sst       nao              0.5122988
## 185      fbar    totaln             -0.5075128
## 245     chel1       nao             -0.5058363

Wykres prezentujący graficznie korelacje między wszystkimi zmiennymi prezentuje się następująco:

corrplot(cor, type="upper")

Można zaobserwować wysoką korelację między dostępnością różnych rodzajów planktonu. Spowodowane jest to prawdopodobnie tym, że pomimo różnic międzygatunkowych każdy rodzaj planktonu potrzebuje mniej więcej zbliżonych warunków do życia i czynniki które wpływają na wzrost lub spadek liczby planktonu jednego rodzaju wpływają w analogiczny sposób na występowanie innych gatunków. Część korelacji jest stosunkowo oczywista i wynika z powiązania między atrybutami. Przykładowo zmienne fbar i cumf dotyczą tej samej wartości (natężenia połowów w regionie) i występuje między nimi bezpośrednia zależność. Atrybutem najmocniej skorelowanym z długością śledzi jest sst oznaczający temperaturę przy powierzchni wody.

Zmiana rozmiaru śledzi w czasie

Bardzo duża liczba pobranych wyników stawia wyzwanie przed czytelnym przedstawieniem zmiany rozmiaru śledzia w czasie. Pierwszy z pokazanych wykresów dokonuje uśrednienia 200 obserwacji występujących obok siebie:

data_len <- select(data, X, length)
n <- 200
data_aggr <- aggregate(data_len,list(rep(1:(nrow(data_len)%/%n+1),each=n,len=nrow(data_len))),mean)[-1]
ggplotly(ggplot(data_aggr, aes(x=X, y=length)) + geom_line())

Drugie podejście korzysta z obecnej w paczce ggplot metody geom_smooth, która wygładza wykres. W tym przypadku jako argument funkcji przekazywany jest cały zbiór:

ggplotly(ggplot(data, aes(x=X, y=length)) + geom_smooth())
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Oba wykresy pozwalają zauważyć, że rozmiar śledzi przez pierwsze około 17 tysięcy pomiarów, a następnie stopniowo spada.

Regresor przewidujący rozmiar śledzia

Aby móc porównywać między sobą różne modele regresji zdefiniowano 5-krotną walidację krzyżową powtórzoną dwukrotnie. Dodatkowo przed uczeniem modelu wartości atrybutów są normalizowane poprzez odjęcie wartości średniej wewnątrz atrybutu oraz podzielenie wyniku przez odchylenie standardowe:

fitControl <- trainControl(method = "repeatedcv", number = 5, repeats = 2, preProcOptions = c('scale', 'center'))

Celem analizy jest zbadanie jaki wpływ na rozmiar śledzi ma środowisko, w którym żyją. Z tego powodu modele regresji trenowane są na danych pozbawionych atrybutów związanych z czasem, czyli “X” oraz “xmonth”:

dataLearning <- select(data, -X, -xmonth)

W ramach analizy skonstruowano modele regresji korzystając z wielu metod: * Regresja liniowa:

modelLm <- train(length ~ ., data = dataLearning, method='lm', trControl=fitControl)
modelLm
## Linear Regression 
## 
## 52582 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times) 
## Summary of sample sizes: 42063, 42066, 42066, 42067, 42066, 42066, ... 
## Resampling results:
## 
##   RMSE      Rsquared   MAE    
##   1.362265  0.3207204  1.08217
## 
## Tuning parameter 'intercept' was held constant at a value of TRUE
modelLars <- train(length ~ ., data = dataLearning, method='lars', trControl=fitControl)
modelLars
## Least Angle Regression 
## 
## 52582 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times) 
## Summary of sample sizes: 42064, 42066, 42065, 42066, 42067, 42066, ... 
## Resampling results across tuning parameters:
## 
##   fraction  RMSE      Rsquared   MAE     
##   0.050     1.576921  0.2046328  1.273636
##   0.525     1.381235  0.3059592  1.100364
##   1.000     1.362400  0.3206643  1.082284
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was fraction = 1.
modelRidge <- train(length ~ ., data = dataLearning, method='ridge', trControl=fitControl)
modelRidge
## Ridge Regression 
## 
## 52582 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times) 
## Summary of sample sizes: 42065, 42066, 42065, 42066, 42066, 42065, ... 
## Resampling results across tuning parameters:
## 
##   lambda  RMSE      Rsquared   MAE     
##   0e+00   1.362488  0.3205524  1.082218
##   1e-04   1.362488  0.3205528  1.082224
##   1e-01   1.375899  0.3074362  1.095358
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was lambda = 1e-04.
modelRf <- train(length ~ ., data = dataLearning, method='rf', ntree=25, importance=T, trControl=fitControl)
modelRf
## Random Forest 
## 
## 52582 samples
##    13 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 2 times) 
## Summary of sample sizes: 42066, 42066, 42066, 42065, 42065, 42066, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE      
##    2    1.188379  0.4831498  0.9407136
##    7    1.194283  0.4780936  0.9447520
##   13    1.196616  0.4761352  0.9466780
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.

Analiza ważności atrybutów

Najlepszą trafność uzyskano wykorzystując regresor Random Forest. Wyliczona ważność atrybutów prezentuje się następująco:

varImp(modelRf)
## rf variable importance
## 
##        Overall
## chel1  100.000
## lcop2   75.653
## sst     61.480
## cfin2   58.180
## recr    45.469
## totaln  41.738
## chel2   26.446
## nao     22.441
## fbar    22.142
## cumf    21.827
## sal     19.187
## lcop1    3.233
## cfin1    0.000

Przedstawione wyniki są znormalizowane, co oznacza, że najbardziej znaczący atrybut posiada wartość 100, a najmniej znaczący wartość 0. Atrybutem, który w największym stopniu odpowiada za zmianę długości śledzia jest atrybut cfin2 oznaczający dostępność planktonu Calanus finmarchicus gat. 2. Trzy kolejne atrybuty w kolejności znaczenia to totaln (łączna liczba ryb złowionych w ramach połowu), sst (temperatura przy powierzchni wody) oraz fbar (natężenie połowów w regionie). Wykresy dla tych atrybutów prezentują się następująco:

plots_cfin2 <- ggplot(data, aes(x=X, y=cfin2)) + geom_smooth()
plots_totaln <- ggplot(data, aes(x=X, y=totaln)) + geom_smooth()
plots_sst <- ggplot(data, aes(x=X, y=sst)) + geom_smooth()
plots_fbar <- ggplot(data, aes(x=X, y=fbar)) + geom_smooth()

grid.arrange(plots_cfin2, plots_totaln,plots_sst, plots_fbar, ncol=2)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Pierwszym wnioskiem, który można wysunąć jest podobieństwo kształtów wykresów dla atrybutów totaln oraz fbar. Atrybut fbar oznacza natężenie połowów totaln łączną liczbę śledzi. Wykresy potwierdzają intuicję, która podpowiada, że zwiększone połowy będą skutkować obniżeniem liczby śledzi zwłowionych w ramach pojedynczego połowu. Zjawisko to można zaobserwować np. podczas pierwszych 10 tysięcy obserwacji, gdy meljąca intensyfikacja połowów skutkowała dynamicznym wzrostem liczby złowionych jednorazowo ryb.

W przypadku ptzebiegu atrybutu cfn2 możemy zaobserwować przebieg zbliżony do przebiegu długości śledzi, lecz z lekkim przesunięciem. Atrybut ten oznacza występowanie planktonu konkretnego gatunku. Możliwym wyjaśnieniem obserwowanych zmian jest fakt, że rosnąca liczba ryb powoduje zmniejszenie ilości dostępnego planktonu, który stanowi dla nich pożywienie. Tłumaczy to początkowy spadek ilości planktonu, gdy liczba śledzi rośnie oraz późniejszy wzrost, gdy liczba ryb maleje. Tajemniczy pozostaje jednak gwałtowny spadek następujący mniej więcej po 35 tysiącach pomiarów. Atrybut sst charakteryzuje dokładnie odwrotny przebieg - początkowy spadek zmienia się w połowie pomiarów w znaczący wzrost. Analiza korelacji pokazała dodatkowo, że między atrybutem sst oraz długością śledzia występuje silna znacząca korelacja ujemna. Pozwala to sądzić, że przyczyną, która w dużym stopniu spowodowała spadek liczby śledzi może być wzrost temperatury przy powierzchni wody.